#### Meta Data ####

# article: "The Ambivalent Effect of Autocratization on Domestic Terrorism" 
# journal: "Studies in Conflict & Terrorism"
# authors: "Lott, Lars; Croissant, Aurel; Trinn. Christoph"
# date: 2023-10-02
# written under "R version 4.3.1 (2023-06-16 ucrt)"


#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# set working directory
getwd()

#### Preliminaries ####

library(plm)
library(MASS)
library(multiwayvcov)
library(ggeffects)
library(PanelMatch)
library(tidyverse)
library(ggpubr)
library(here)
library(haven)
library(imputeTS)
library(modelsummary)


################################################################################
########################## Supplementary Appendix B ############################
################################################################################

#### Load Dataset ####

vdem <- readRDS("data/vdem12.rds") 

vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(lag_v2x_polyarchy = dplyr::lag(v2x_polyarchy, 1), 
         delta_v2x_polyarchy = v2x_polyarchy - lag_v2x_polyarchy) %>%
  mutate(autocratization_democracy = ifelse(delta_v2x_polyarchy <= -0.1 & dplyr::lag(v2x_regime) >=2, 1, 0), 
         autocratization_autocracy = ifelse(delta_v2x_polyarchy <= -0.1 & dplyr::lag(v2x_regime) < 2, 1, 0), 
         democratization_democracy = ifelse(delta_v2x_polyarchy >= 0.1 & dplyr::lag(v2x_regime) >=2, 1, 0), 
         democratization_autocracy = ifelse(delta_v2x_polyarchy >= 0.1 & dplyr::lag(v2x_regime) <2, 1, 0))

vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(autocratization_treatment = ifelse(autocratization_democracy == 1 | autocratization_autocracy == 1, 1, 0))

table(vdem$autocratization_democracy)
table(vdem$autocratization_autocracy)
table(vdem$democratization_democracy)
table(vdem$democratization_autocracy)
table(vdem$autocratization_treatment)

example_data <- vdem %>% 
  filter(year >= 1950) %>%
  dplyr::select(country_name, year, v2x_polyarchy, lag_v2x_polyarchy, delta_v2x_polyarchy, v2x_regime,  
                autocratization_democracy, autocratization_autocracy, democratization_democracy, 
                democratization_autocracy, autocratization_treatment)

# Further Defining Treatments #

vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(democratic_breakdown_delta = ifelse(autocratization_democracy == 1 & (dplyr::lag(v2x_regime) >=2 | dplyr::lag(v2x_regime, 2) >=2) & v2x_regime <2, 1, 0), 
         democratic_breakdown_delta = ifelse(autocratization_autocracy == 1 & (dplyr::lag(v2x_regime) >=2 | dplyr::lag(v2x_regime, 2) >=2) & v2x_regime <2, 1, democratic_breakdown_delta), 
         democratic_breakdown_delta = ifelse(dplyr::lag(democratic_breakdown_delta) == 1 & autocratization_autocracy == 1, 1, democratic_breakdown_delta), 
         democratic_breakdown_delta = ifelse(dplyr::lead(democratic_breakdown_delta) == 1 & autocratization_democracy == 1, 1, democratic_breakdown_delta), 
         democratic_breakdown_delta = ifelse((dplyr::lead(v2x_regime) <2 | dplyr::lead(v2x_regime, 2) <2) & autocratization_democracy == 1, 1, democratic_breakdown_delta), 
         autocratization_autocracy = ifelse(democratic_breakdown_delta == 1, 0, autocratization_autocracy), 
         autocratization_democracy = ifelse(democratic_breakdown_delta == 1, 0, autocratization_democracy))

table(vdem$democratic_breakdown_delta) # 72 country-years
table(vdem$autocratization_autocracy) # 94 country-years
table(vdem$autocratization_democracy) # 4 country-years 


vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(democratic_transition_delta = ifelse(democratization_autocracy == 1 & (v2x_regime >=2 | dplyr::lead(v2x_regime) >= 2), 1, 0), 
         democratic_transition_delta = ifelse(democratization_democracy == 1 & (dplyr::lag(v2x_regime) <2 | dplyr::lag(v2x_regime, 2) <2 ) & v2x_regime>= 2, 1, democratic_transition_delta), 
         democratic_transition_delta = ifelse(dplyr::lead(democratic_transition_delta) == 1 & democratization_autocracy == 1, 1, democratic_transition_delta),
         democratic_transition_delta = ifelse(dplyr::lag(democratic_transition_delta) == 1 & democratization_democracy == 1, 1, democratic_transition_delta),
         democratization_democracy = ifelse(democratic_transition_delta == 1, 0, democratization_democracy), 
         democratization_autocracy = ifelse(democratic_transition_delta == 1, 0, democratization_autocracy))

example_data <- vdem %>% 
  filter(year >= 1950) %>%
  dplyr::select(country_name, year, v2x_polyarchy, lag_v2x_polyarchy, delta_v2x_polyarchy, v2x_regime,  
                autocratization_democracy, democratic_breakdown_delta, autocratization_autocracy, democratization_democracy, 
                democratic_transition_delta, democratization_autocracy)

table(vdem$democratic_transition_delta) # 172 country-years
table(vdem$democratization_democracy) # 20 country-years
table(vdem$democratization_autocracy) # 164 country-years 



vdem <- vdem %>%
  mutate(log_count_attacks_dom = log(count_attacks_dom +0.01),
         log_count_attacks_int = log(count_attacks_int +0.01), 
         log_count_attacks = log(count_attacks +0.01),
         log_sum_nkill_dom = log(sum_nkill_dom +0.01),
         log_sum_nkill_int = log(sum_nkill_int +0.01),
         log_sum_nkill = log(sum_nkill +0.01)) %>%
  filter(year >= 1966)


vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(log_count_attacks_dom = na_interpolation(log_count_attacks_dom, option = "spline",  maxgap = 1), 
         log_sum_nkill_dom = na_interpolation(log_sum_nkill_dom, option = "spline",  maxgap = 1), 
         count_attacks_dom = na_interpolation(count_attacks_dom, option = "spline",  maxgap = 1), 
         sum_nkill_dom = na_interpolation(sum_nkill_dom, option = "spline",  maxgap = 1), 
         log_count_attacks_dom = ifelse(log_count_attacks_dom<0, 0, log_count_attacks_dom),
         log_sum_nkill_dom = ifelse(log_sum_nkill_dom<0, 0, log_sum_nkill_dom),
         count_attacks_dom = ifelse(count_attacks_dom<0, 0, count_attacks_dom),
         sum_nkill_dom = ifelse(sum_nkill_dom<0, 0, sum_nkill_dom))


vdem <- vdem %>%
  dplyr::select(country_id, country_name, year, aut_ep, aut_ep_id, 
                aut_ep_id,  log_count_attacks_dom, log_count_attacks_int, log_count_attacks, 
                log_sum_nkill_dom, log_sum_nkill_int, log_sum_nkill, count_attacks_dom, count_attacks_int, count_attacks, 
                sum_nkill_dom, sum_nkill_int, sum_nkill,e_gdppc, e_pop,
                v2x_polyarchy, milper_extrapol, eprratio_extrapol, 
                civil_war_ucdp, v2x_regime, aut_ep_outcome, aut_ep_outcome_agg, 
                v2caviol, dem_ep, dem_ep_id, dem_ep_outcome_agg, dem_ep_outcome, e_regionpol, e_regionpol_6C, 
                autocratization_democracy, democratic_breakdown_delta, autocratization_autocracy, democratization_democracy, 
                democratic_transition_delta, democratization_autocracy, lag_v2x_polyarchy, delta_v2x_polyarchy, 
                autocratization_treatment, ma_count_attacks_dom, ma_sum_nkill_dom, democracy_stock, 
                v2x_veracc, v2x_diagacc, v2x_horacc, e_pelifeex, e_peaveduc, e_area) %>%
  filter(year <=2020) %>%
  group_by(country_id) %>%
  fill(e_area)

vdem <- vdem %>%
  drop_na(autocratization_democracy) %>%
  mutate(e_gdppc_sqaured = e_gdppc^2, 
         e_area_log = log10(e_area))

summary(vdem$e_area_log)
summary(vdem$e_gdppc_sqaured)



#------------------------------------------------------------------------------#
#------------------------------ Figure B1        ------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)


vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$democratic_breakdown_delta)
table(vdem$democratic_breakdown_delta)


#### 1. Effect of democratic_breakdown_delta on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "democratic_breakdown_delta", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_5, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "democratic_breakdown_delta", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_5, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "democratic_breakdown_delta", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_5, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "democratic_breakdown_delta", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets) # 32 democratic breakdown

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "democratic_breakdown_delta", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_B1/01_Figure_Frequency_Distribution.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for democratic_breakdown_delta", rot = 90))
ggsave("results/Figure_B1/01_Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_B1/01_Figure_Improved_Covariate_Balance.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "democratic_breakdown_delta", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "democratic_breakdown_delta", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))


saveRDS(M1_pred, "results/Figure_B1/data_M1_pred_dem_breakdown.rds")
saveRDS(M2_pred, "results/Figure_B1/data_M2_pred_dem_breakdown.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

saveRDS(M1_pred_ci09, "results/Figure_B1/data_M1_pred_dem_breakdown_ci09.rds")
saveRDS(M2_pred_ci09, "results/Figure_B1/data_M2_pred_dem_breakdown_ci09.rds")


#### Regressed Autocracy Effect ####

data_figure_1a <- readRDS("results/Figure_B1/data_M1_pred_dem_breakdown.rds") 
data_figure_1b <- readRDS("results/Figure_B1/data_M2_pred_dem_breakdown.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_B1/data_M1_pred_dem_breakdown_ci09.rds") 
data_figure_1b_09 <- readRDS("results/Figure_B1/data_M2_pred_dem_breakdown_ci09.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = time)) +
  geom_linerange(aes(ymin=`2.5%`, ymax = `97.5%`), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=`5%`, ymax = `95%`), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(0,10,1)) +
  labs(x = "Year around Democratic Breakdown/Recession", 
       y = "Estimated effect of Democratic Breakdown/Recession", 
       title = "Propensity Score Weighting: Terrorist Attacks")

f1a 

ggsave("outputs/Figure_B1.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_B1.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)


#------------------------------------------------------------------------------#
#-------------------------------- Figure B2   ---------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)

vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$autocratization_autocracy)


#### 1. Effect of Autocratization on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "autocratization_autocracy", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_5, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "autocratization_autocracy", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_5, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "autocratization_autocracy", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_5, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "autocratization_autocracy", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets)

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "autocratization_autocracy", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_B2/01_Figure_Frequency_Distribution.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for Autocratization", rot = 90))
ggsave("results/Figure_B2/01_Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_B2/01_Figure_Improved_Covariate_Balance.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_autocracy", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_autocracy", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))


saveRDS(M1_pred, "results/Figure_B2/data_M1_pred_dem_breakdown.rds")
saveRDS(M2_pred, "results/Figure_B2/data_M2_pred_dem_breakdown.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

saveRDS(M1_pred_ci09, "results/Figure_B2/data_M1_pred_dem_breakdown_ci09.rds")
saveRDS(M2_pred_ci09, "results/Figure_B2/data_M2_pred_dem_breakdown_ci09.rds")


#### Regressed Autocracy Effect ####

data_figure_1a <- readRDS("results/Figure_B2/data_M1_pred_dem_breakdown.rds") 
data_figure_1b <- readRDS("results/Figure_B2/data_M2_pred_dem_breakdown.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_B2/data_M1_pred_dem_breakdown_ci09.rds") 
data_figure_1b_09 <- readRDS("results/Figure_B2/data_M2_pred_dem_breakdown_ci09.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = time)) +
  geom_linerange(aes(ymin=`2.5%`, ymax = `97.5%`), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=`5%`, ymax = `95%`), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(0,10,1)) +
  labs(x = "Year around Autocratic Consolidation", 
       y = "Estimated effect of Autocratic Consolidation", 
       title = "Propensity Score Weighting: Terrorist Attacks")
f1a

ggsave("outputs/Figure_B2.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_B2.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)


#------------------------------------------------------------------------------#
#--------------------------------- Figure B3 ----------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)

vdem <- vdem %>%
  dplyr::select(-c(country_name, aut_ep_id))

vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$autocratization_treatment)
table(vdem_subset$autocratization_treatment)

table(vdem_subset$democratization_autocracy)
table(vdem_subset$democratic_transition_delta)


min(vdem$year)
max(vdem$year)

#### 1. Effect of democratic_breakdown_delta on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "autocratization_treatment", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_5, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "autocratization_treatment", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_5, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "autocratization_treatment", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_5, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets) # 63 treatments

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "autocratization_treatment", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_B3/01_Figure_Frequency_Distribution.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for democratic_breakdown_delta", rot = 90))
ggsave("results/Figure_B3/01_Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom",  "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom",  "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom",  "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_B3/01_Figure_Improved_Covariate_Balance.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                          I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))


saveRDS(M1_pred, "results/Figure_B3/data_M1_pred_dem_breakdown.rds")
saveRDS(M2_pred, "results/Figure_B3/data_M2_pred_dem_breakdown.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

saveRDS(M1_pred_ci09, "results/Figure_B3/data_M1_pred_dem_breakdown_ci09.rds")
saveRDS(M2_pred_ci09, "results/Figure_B3/data_M2_pred_dem_breakdown_ci09.rds")


#### Autocratization Effect ####

data_figure_1a <- readRDS("results/Figure_B3/data_M1_pred_dem_breakdown.rds") 
data_figure_1b <- readRDS("results/Figure_B3/data_M2_pred_dem_breakdown.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_B3/data_M1_pred_dem_breakdown_ci09.rds") 
data_figure_1b_09 <- readRDS("results/Figure_B3/data_M2_pred_dem_breakdown_ci09.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = time)) +
  geom_linerange(aes(ymin=`2.5%`, ymax = `97.5%`), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=`5%`, ymax = `95%`), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(0,10,1)) +
  labs(x = "Year around Autocraitzation", 
       y = "Estimated effect of Autocraitzation", 
       title = "Propensity Score Weighting: Terrorist Attacks")

f1a
ggsave("outputs/Figure_B3.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_B3.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)

################################################################################
########################## Supplementary Appendix C ############################
################################################################################


#### Load Dataset ####

vdem <- readRDS("data/vdem12.rds") 

#### Define Treatment ####

vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(lag_v2x_polyarchy = dplyr::lag(v2x_polyarchy, 1), 
         delta_v2x_polyarchy = v2x_polyarchy - lag_v2x_polyarchy) %>%
  mutate(autocratization_democracy = ifelse(delta_v2x_polyarchy <= -0.05 & dplyr::lag(v2x_regime) >=2, 1, 0), 
         autocratization_autocracy = ifelse(delta_v2x_polyarchy <= -0.05 & dplyr::lag(v2x_regime) < 2, 1, 0), 
         democratization_democracy = ifelse(delta_v2x_polyarchy >= 0.05 & dplyr::lag(v2x_regime) >=2, 1, 0), 
         democratization_autocracy = ifelse(delta_v2x_polyarchy >= 0.05 & dplyr::lag(v2x_regime) <2, 1, 0))

vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(autocratization_treatment = ifelse(autocratization_democracy == 1 | autocratization_autocracy == 1, 1, 0))

table(vdem$autocratization_democracy)
table(vdem$autocratization_autocracy)
table(vdem$democratization_democracy)
table(vdem$democratization_autocracy)
table(vdem$autocratization_treatment)

example_data <- vdem %>% 
    filter(year >= 1950) %>%
    dplyr::select(country_name, year, v2x_polyarchy, lag_v2x_polyarchy, delta_v2x_polyarchy, v2x_regime,  
                  autocratization_democracy, autocratization_autocracy, democratization_democracy, 
                  democratization_autocracy, autocratization_treatment)

# Further Defining Treatments #

vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(democratic_breakdown_delta = ifelse(autocratization_democracy == 1 & (dplyr::lag(v2x_regime) >=2 | dplyr::lag(v2x_regime, 2) >=2) & v2x_regime <2, 1, 0), 
         democratic_breakdown_delta = ifelse(autocratization_autocracy == 1 & (dplyr::lag(v2x_regime) >=2 | dplyr::lag(v2x_regime, 2) >=2) & v2x_regime <2, 1, democratic_breakdown_delta), 
         democratic_breakdown_delta = ifelse(dplyr::lag(democratic_breakdown_delta) == 1 & autocratization_autocracy == 1, 1, democratic_breakdown_delta), 
         democratic_breakdown_delta = ifelse(dplyr::lead(democratic_breakdown_delta) == 1 & autocratization_democracy == 1, 1, democratic_breakdown_delta), 
         democratic_breakdown_delta = ifelse((dplyr::lead(v2x_regime) <2 | dplyr::lead(v2x_regime, 2) <2) & autocratization_democracy == 1, 1, democratic_breakdown_delta), 
         autocratization_autocracy = ifelse(democratic_breakdown_delta == 1, 0, autocratization_autocracy), 
         autocratization_democracy = ifelse(democratic_breakdown_delta == 1, 0, autocratization_democracy))

table(vdem$democratic_breakdown_delta) # 65 country-years
table(vdem$autocratization_autocracy) # 100 country-years
table(vdem$autocratization_democracy) # 8 country-years 


vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(democratic_transition_delta = ifelse(democratization_autocracy == 1 & (v2x_regime >=2 | dplyr::lead(v2x_regime) >= 2), 1, 0), 
         democratic_transition_delta = ifelse(democratization_democracy == 1 & (dplyr::lag(v2x_regime) <2 | dplyr::lag(v2x_regime, 2) <2 ) & v2x_regime>= 2, 1, democratic_transition_delta), 
         democratic_transition_delta = ifelse(dplyr::lead(democratic_transition_delta) == 1 & democratization_autocracy == 1, 1, democratic_transition_delta),
         democratic_transition_delta = ifelse(dplyr::lag(democratic_transition_delta) == 1 & democratization_democracy == 1, 1, democratic_transition_delta),
         democratization_democracy = ifelse(democratic_transition_delta == 1, 0, democratization_democracy), 
         democratization_autocracy = ifelse(democratic_transition_delta == 1, 0, democratization_autocracy))

example_data <- vdem %>% 
  filter(year >= 1950) %>%
  dplyr::select(country_name, year, v2x_polyarchy, lag_v2x_polyarchy, delta_v2x_polyarchy, v2x_regime,  
                autocratization_democracy, democratic_breakdown_delta, autocratization_autocracy, democratization_democracy, 
                democratic_transition_delta, democratization_autocracy)

table(vdem$democratic_transition_delta) # 145 country-years
table(vdem$democratization_democracy) # 38 country-years
table(vdem$democratization_autocracy) # 173 country-years 



vdem <- vdem %>%
  mutate(log_count_attacks_dom = log(count_attacks_dom +0.01),
         log_count_attacks_int = log(count_attacks_int +0.01), 
         log_count_attacks = log(count_attacks +0.01),
         log_sum_nkill_dom = log(sum_nkill_dom +0.01),
         log_sum_nkill_int = log(sum_nkill_int +0.01),
         log_sum_nkill = log(sum_nkill +0.01)) %>%
  filter(year >= 1966)


vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(log_count_attacks_dom = na_interpolation(log_count_attacks_dom, option = "spline",  maxgap = 1), 
         log_sum_nkill_dom = na_interpolation(log_sum_nkill_dom, option = "spline",  maxgap = 1), 
         count_attacks_dom = na_interpolation(count_attacks_dom, option = "spline",  maxgap = 1), 
         sum_nkill_dom = na_interpolation(sum_nkill_dom, option = "spline",  maxgap = 1), 
         log_count_attacks_dom = ifelse(log_count_attacks_dom<0, 0, log_count_attacks_dom),
         log_sum_nkill_dom = ifelse(log_sum_nkill_dom<0, 0, log_sum_nkill_dom),
         count_attacks_dom = ifelse(count_attacks_dom<0, 0, count_attacks_dom),
         sum_nkill_dom = ifelse(sum_nkill_dom<0, 0, sum_nkill_dom))


vdem <- vdem %>%
  dplyr::select(country_id, country_name, year, aut_ep, aut_ep_id, 
                aut_ep_id,  log_count_attacks_dom, log_count_attacks_int, log_count_attacks, 
                log_sum_nkill_dom, log_sum_nkill_int, log_sum_nkill, count_attacks_dom, count_attacks_int, count_attacks, 
                sum_nkill_dom, sum_nkill_int, sum_nkill,e_gdppc, e_pop,
                v2x_polyarchy, milper_extrapol, eprratio_extrapol, 
                civil_war_ucdp, v2x_regime, aut_ep_outcome, aut_ep_outcome_agg, 
                v2caviol, dem_ep, dem_ep_id, dem_ep_outcome_agg, dem_ep_outcome, e_regionpol, e_regionpol_6C, 
                autocratization_democracy, democratic_breakdown_delta, autocratization_autocracy, democratization_democracy, 
                democratic_transition_delta, democratization_autocracy, lag_v2x_polyarchy, delta_v2x_polyarchy, 
                autocratization_treatment, ma_count_attacks_dom, ma_sum_nkill_dom, democracy_stock, 
                v2x_veracc, v2x_diagacc, v2x_horacc, e_pelifeex, e_peaveduc, e_area) %>%
  filter(year <=2020) %>%
  group_by(country_id) %>%
  fill(e_area)

vdem <- vdem %>%
  drop_na(autocratization_democracy) %>%
  mutate(e_gdppc_sqaured = e_gdppc^2, 
         e_area_log = log10(e_area))

summary(vdem$e_area_log)
summary(vdem$e_gdppc_sqaured)


#------------------------------------------------------------------------------#
#------------------------- Deskriptive Tables for alle QoI -------------------#
#------------------------------------------------------------------------------#

vdem_subset <- vdem %>%
  filter(year >= 1970) 

# Democratic Breakdown #
list_countries_dem_breakdown <- vdem_subset %>%
  group_by(country_name) %>%
  filter(democratic_breakdown_delta == 1) %>%
  dplyr::select(country_name, year, democratic_breakdown_delta, v2x_polyarchy, lag_v2x_polyarchy) %>%
  mutate(year = as.factor(year), 
         democratic_breakdown_delta = as.factor(democratic_breakdown_delta))

datasummary_df(list_countries_dem_breakdown)

# Regressed Autocracy #
list_countries_regressed_autocracy <- vdem_subset %>%
  group_by(country_name) %>%
  filter(autocratization_autocracy == 1) %>%
  dplyr::select(country_name, year, autocratization_autocracy, v2x_polyarchy, lag_v2x_polyarchy) %>%
  mutate(year = as.factor(year), 
         autocratization_autocracy = as.factor(autocratization_autocracy))

datasummary_df(list_countries_regressed_autocracy)

# Democratic Regression #
list_countries_democratic_regression <- vdem_subset %>%
  group_by(country_name) %>%
  filter(autocratization_democracy == 1) %>%
  dplyr::select(country_name, year, autocratization_democracy, v2x_polyarchy, lag_v2x_polyarchy) %>%
  mutate(year = as.factor(year), 
         autocratization_democracy = as.factor(autocratization_democracy))

datasummary_df(list_countries_democratic_regression)

# Democratic Transition #
list_countries_democratic_transition <- vdem_subset %>%
  group_by(country_name) %>%
  filter(democratic_transition_delta == 1) %>%
  dplyr::select(country_name, year, democratic_transition_delta, v2x_polyarchy, lag_v2x_polyarchy) %>%
  mutate(year = as.factor(year), 
         democratic_transition_delta = as.factor(democratic_transition_delta))

datasummary_df(list_countries_democratic_transition )

# Liberalized Autocracy #
list_countries_liberalized_autocracy <- vdem_subset %>%
  group_by(country_name) %>%
  filter(democratization_autocracy == 1) %>%
  dplyr::select(country_name, year, democratization_autocracy, v2x_polyarchy, lag_v2x_polyarchy) %>%
  mutate(year = as.factor(year), 
         democratization_autocracy = as.factor(democratization_autocracy))

datasummary_df(list_countries_liberalized_autocracy)

# Liberalized Democracy #
list_countries_liberalized_democracy <- vdem_subset %>%
  group_by(country_name) %>%
  filter(democratization_democracy == 1) %>%
  dplyr::select(country_name, year, democratization_democracy, v2x_polyarchy, lag_v2x_polyarchy) %>%
  mutate(year = as.factor(year), 
         democratization_democracy = as.factor(democratization_democracy))

datasummary_df(list_countries_liberalized_democracy)


#------------------------------------------------------------------------------#
#--------------------------------- Figure C1 ----------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)


vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$democratic_breakdown_delta)
table(vdem$democratic_breakdown_delta)


#### 1. Effect of democratic_breakdown_delta on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "democratic_breakdown_delta", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_5, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "democratic_breakdown_delta", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_5, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "democratic_breakdown_delta", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_5, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "democratic_breakdown_delta", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets) # 44 treatments

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "democratic_breakdown_delta", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_C1/Figure_Frequency_Distribution.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for democratic_breakdown_delta", rot = 90))
ggsave("results/Figure_C1/Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_C1/Figure_Improved_Covariate_Balance.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "democratic_breakdown_delta", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "democratic_breakdown_delta", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####
 
#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))


saveRDS(M1_pred, "results/Figure_C1/data_M1_pred_dem_breakdown.rds")
saveRDS(M2_pred, "results/Figure_C1/data_M2_pred_dem_breakdown.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

saveRDS(M1_pred_ci09, "results/Figure_C1/data_M1_pred_dem_breakdown_ci09.rds")
saveRDS(M2_pred_ci09, "results/Figure_C1/data_M2_pred_dem_breakdown_ci09.rds")


#### Democratic Breakdown Effect ####

data_figure_1a <- readRDS("results/Figure_C1/data_M1_pred_dem_breakdown.rds") 
data_figure_1b <- readRDS("results/Figure_C1/data_M2_pred_dem_breakdown.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_C1/data_M1_pred_dem_breakdown_ci09.rds") 
data_figure_1b_09 <- readRDS("results/Figure_C1/data_M2_pred_dem_breakdown_ci09.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = time)) +
  geom_linerange(aes(ymin=`2.5%`, ymax = `97.5%`), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=`5%`, ymax = `95%`), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(0,10,1)) +
  labs(x = "Year around Democratic Breakdown", 
       y = "Estimated effect of Democratic Breakdown", 
       title = "Propensity Score Weighting: Terrorist Attacks")


f1a 

ggsave("outputs/Figure_C1.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_C1.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)

#------------------------------------------------------------------------------#
#-------------------------------- Figure C2            ------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)

vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$autocratization_autocracy)


#### 1. Effect of Autocratization on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "autocratization_democracy", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_5, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "autocratization_democracy", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_5, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "autocratization_democracy", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_5, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "autocratization_democracy", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets) # 11 treatments

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "autocratization_democracy", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_C2/Figure_Frequency_Distribution.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for Autocratization", rot = 90))
ggsave("results/Figure_C2/Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_C2/Figure_Improved_Covariate_Balance.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_democracy", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_democracy", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))


saveRDS(M1_pred, "results/Figure_C2/data_M1_pred_dem_breakdown.rds")
saveRDS(M2_pred, "results/Figure_C2/data_M2_pred_dem_breakdown.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

saveRDS(M1_pred_ci09, "results/Figure_C2/data_M1_pred_dem_breakdown_ci09.rds")
saveRDS(M2_pred_ci09, "results/Figure_C2/data_M2_pred_dem_breakdown_ci09.rds")


#### Regressed Autocracy Effect ####

data_figure_1a <- readRDS("results/Figure_C2/data_M1_pred_dem_breakdown.rds") 
data_figure_1b <- readRDS("results/Figure_C2/data_M2_pred_dem_breakdown.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_C2/data_M1_pred_dem_breakdown_ci09.rds") 
data_figure_1b_09 <- readRDS("results/Figure_C2/data_M2_pred_dem_breakdown_ci09.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = time)) +
  geom_linerange(aes(ymin=`2.5%`, ymax = `97.5%`), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=`5%`, ymax = `95%`), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(0,10,1)) +
  labs(x = "Year around autocratization", 
       y = "Estimated effect of Autocratization", 
       title = "Propensity Score Weighting: Terrorist Attacks")

f1a 
ggsave("outputs/Figure_C2.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_C2.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)


#------------------------------------------------------------------------------#
#--------------------------------- Figure C3     ------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)

vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$autocratization_autocracy)


#### 1. Effect of Autocratization on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "autocratization_autocracy", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_5, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "autocratization_autocracy", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_5, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "autocratization_autocracy", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_5, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "autocratization_autocracy", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets)

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "autocratization_autocracy", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_C3/Figure_Frequency_Distribution.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for Autocratization", rot = 90))
ggsave("results/Figure_C3/Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_C3/Figure_Improved_Covariate_Balance.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_autocracy", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)),  
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_autocracy", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))


saveRDS(M1_pred, "results/Figure_C3/data_M1_pred_dem_breakdown.rds")
saveRDS(M2_pred, "results/Figure_C3/data_M2_pred_dem_breakdown.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

saveRDS(M1_pred_ci09, "results/Figure_C3/data_M1_pred_dem_breakdown_ci09.rds")
saveRDS(M2_pred_ci09, "results/Figure_C3/data_M2_pred_dem_breakdown_ci09.rds")


#### Regressed Autocracy Effect ####

data_figure_1a <- readRDS("results/Figure_C3/data_M1_pred_dem_breakdown.rds") 
data_figure_1b <- readRDS("results/Figure_C3/data_M2_pred_dem_breakdown.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_C3/data_M1_pred_dem_breakdown_ci09.rds") 
data_figure_1b_09 <- readRDS("results/Figure_C3/data_M2_pred_dem_breakdown_ci09.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = time)) +
  geom_linerange(aes(ymin=`2.5%`, ymax = `97.5%`), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=`5%`, ymax = `95%`), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(0,10,1)) +
  labs(x = "Year around Autocratic Consolidation", 
       y = "Estimated effect of Autocratic Consolidation", 
       title = "Propensity Score Weighting: Terrorist Attacks")

f1a
ggsave("outputs/Figure_C3.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_C3.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)


#------------------------------------------------------------------------------#
#--------------------------------- Figure C4 ----------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_5 <- c(0:5)

vdem <- vdem %>%
  dplyr::select(-c(country_name, aut_ep_id))

vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$autocratization_treatment)
table(vdem$autocratization_treatment)

vdem %>%
  filter(year >=1970) %>%
  distinct(country_id)


#### 1. Effect of democratic_breakdown_delta on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "autocratization_treatment", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_5, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "autocratization_treatment", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_5, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "autocratization_treatment", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_5, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_5, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets) # 139 treatments

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "autocratization_treatment", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_C4/Figure_Frequency_Distribution.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for democratic_breakdown_delta", rot = 90))
ggsave("results/Figure_C4/Figure_Covariate_Balance_Autocratization.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_C4/Figure_Improved_Covariate_Balance.pdf', width = 10, height = 5)
dev.off()

##### Autocratizaiton as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                          I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                          I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_5, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5))


saveRDS(M1_pred, "results/Figure_C4/data_M1_pred_dem_breakdown.rds")
saveRDS(M2_pred, "results/Figure_C4/data_M2_pred_dem_breakdown.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5))

saveRDS(M1_pred_ci09, "results/Figure_C4/data_M1_pred_dem_breakdown_ci09.rds")
saveRDS(M2_pred_ci09, "results/Figure_C4/data_M2_pred_dem_breakdown_ci09.rds")


#### Autocratization Effect ####

data_figure_1a <- readRDS("results/Figure_C4/data_M1_pred_dem_breakdown.rds") 
data_figure_1b <- readRDS("results/Figure_C4/data_M2_pred_dem_breakdown.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_C4/data_M1_pred_dem_breakdown_ci09.rds") 
data_figure_1b_09 <- readRDS("results/Figure_C4/data_M2_pred_dem_breakdown_ci09.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = time)) +
  geom_linerange(aes(ymin=`2.5%`, ymax = `97.5%`), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=`5%`, ymax = `95%`), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(0,10,1)) +
  labs(x = "Year around Autocraitzation", 
       y = "Estimated effect of Autocraitzation", 
       title = "Propensity Score Weighting: Terrorist Attacks")

f1a
ggsave("outputs/Figure_C4.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_C4.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)


